Raw data and code from https://drive.google.com/drive/folders/1Jsv34JjNo22NOCd26iBHMP80EPG0xmQE, linked from Casey Handmer's blog post Solar and batteries for generic use cases

In [1]:
from datetime import datetime
from pathlib import Path
from random import uniform

import altair as alt
from numba import njit
import pandas as pd
from vega_datasets import data

alt.data_transformers.enable("vegafusion")

folder = Path('./texas')
folder.exists()
Out[1]:
True
In [2]:
actual = [fp for fp in folder.iterdir() if fp.stem.startswith('Actual')]
actual
Out[2]:
[PosixPath('texas/Actual_36.05_-102.95_2006_UPV_17MW_5_Min.csv'),
 PosixPath('texas/Actual_34.05_-102.05_2006_UPV_151MW_5_Min.csv'),
 PosixPath('texas/Actual_34.95_-101.85_2006_DPV_27MW_5_Min.csv'),
 PosixPath('texas/Actual_33.45_-94.45_2006_DPV_27MW_5_Min.csv'),
 PosixPath('texas/Actual_33.45_-94.35_2006_DPV_27MW_5_Min.csv'),
 PosixPath('texas/Actual_35.85_-102.85_2006_UPV_17MW_5_Min.csv'),
 PosixPath('texas/Actual_35.05_-101.85_2006_DPV_27MW_5_Min.csv'),
 PosixPath('texas/Actual_33.15_-102.25_2006_UPV_50MW_5_Min.csv'),
 PosixPath('texas/Actual_35.35_-101.95_2006_UPV_67MW_5_Min.csv'),
 PosixPath('texas/Actual_32.55_-102.75_2006_UPV_118MW_5_Min.csv'),
 PosixPath('texas/Actual_32.45_-94.85_2006_DPV_35MW_5_Min.csv'),
 PosixPath('texas/Actual_33.15_-94.65_2006_UPV_14MW_5_Min.csv'),
 PosixPath('texas/Actual_34.55_-102.55_2006_UPV_168MW_5_Min.csv'),
 PosixPath('texas/Actual_33.55_-101.85_2006_DPV_37MW_5_Min.csv'),
 PosixPath('texas/Actual_32.35_-94.25_2006_UPV_122MW_5_Min.csv'),
 PosixPath('texas/Actual_33.65_-101.85_2006_DPV_37MW_5_Min.csv'),
 PosixPath('texas/Actual_33.75_-101.95_2006_UPV_101MW_5_Min.csv'),
 PosixPath('texas/Actual_32.55_-94.85_2006_DPV_35MW_5_Min.csv'),
 PosixPath('texas/Actual_32.65_-102.85_2006_UPV_118MW_5_Min.csv'),
 PosixPath('texas/Actual_33.65_-101.75_2006_DPV_37MW_5_Min.csv'),
 PosixPath('texas/Actual_33.05_-94.25_2006_UPV_54MW_5_Min.csv'),
 PosixPath('texas/Actual_32.65_-94.45_2006_UPV_54MW_5_Min.csv'),
 PosixPath('texas/Actual_33.55_-94.45_2006_DPV_27MW_5_Min.csv'),
 PosixPath('texas/Actual_33.75_-102.95_2006_UPV_50MW_5_Min.csv'),
 PosixPath('texas/Actual_32.95_-102.95_2006_UPV_17MW_5_Min.csv'),
 PosixPath('texas/Actual_34.95_-102.95_2006_UPV_67MW_5_Min.csv'),
 PosixPath('texas/Actual_34.15_-101.65_2006_UPV_151MW_5_Min.csv'),
 PosixPath('texas/Actual_32.65_-102.95_2006_UPV_118MW_5_Min.csv'),
 PosixPath('texas/Actual_33.75_-101.85_2006_DPV_37MW_5_Min.csv'),
 PosixPath('texas/Actual_32.75_-102.45_2006_UPV_17MW_5_Min.csv'),
 PosixPath('texas/Actual_32.85_-103.05_2006_UPV_84MW_5_Min.csv'),
 PosixPath('texas/Actual_35.35_-101.75_2006_DPV_28MW_5_Min.csv'),
 PosixPath('texas/Actual_33.45_-101.75_2006_UPV_84MW_5_Min.csv'),
 PosixPath('texas/Actual_35.15_-101.65_2006_UPV_101MW_5_Min.csv'),
 PosixPath('texas/Actual_33.45_-102.95_2006_UPV_151MW_5_Min.csv'),
 PosixPath('texas/Actual_35.75_-102.95_2006_UPV_151MW_5_Min.csv'),
 PosixPath('texas/Actual_33.35_-94.35_2006_UPV_95MW_5_Min.csv'),
 PosixPath('texas/Actual_32.05_-94.15_2006_UPV_95MW_5_Min.csv'),
 PosixPath('texas/Actual_32.05_-94.35_2006_UPV_27MW_5_Min.csv'),
 PosixPath('texas/Actual_35.35_-101.85_2006_DPV_28MW_5_Min.csv'),
 PosixPath('texas/Actual_33.65_-102.45_2006_UPV_168MW_5_Min.csv'),
 PosixPath('texas/Actual_33.65_-101.95_2006_DPV_37MW_5_Min.csv'),
 PosixPath('texas/Actual_32.95_-94.95_2006_UPV_27MW_5_Min.csv'),
 PosixPath('texas/Actual_32.45_-94.75_2006_DPV_35MW_5_Min.csv'),
 PosixPath('texas/Actual_36.25_-102.95_2006_UPV_84MW_5_Min.csv'),
 PosixPath('texas/Actual_34.35_-102.95_2006_UPV_50MW_5_Min.csv'),
 PosixPath('texas/Actual_32.55_-102.25_2006_UPV_118MW_5_Min.csv'),
 PosixPath('texas/Actual_33.55_-94.65_2006_UPV_136MW_5_Min.csv'),
 PosixPath('texas/Actual_33.25_-102.65_2006_UPV_17MW_5_Min.csv'),
 PosixPath('texas/Actual_35.35_-101.95_2006_UPV_17MW_5_Min.csv'),
 PosixPath('texas/Actual_35.45_-101.85_2006_DPV_28MW_5_Min.csv'),
 PosixPath('texas/Actual_32.55_-103.05_2006_UPV_118MW_5_Min.csv'),
 PosixPath('texas/Actual_33.65_-102.05_2006_UPV_34MW_5_Min.csv'),
 PosixPath('texas/Actual_34.65_-102.85_2006_UPV_34MW_5_Min.csv'),
 PosixPath('texas/Actual_34.95_-101.75_2006_DPV_27MW_5_Min.csv'),
 PosixPath('texas/Actual_33.05_-102.95_2006_UPV_67MW_5_Min.csv')]
In [3]:
states = alt.topo_feature(data.us_10m.url, feature='states')

usa_states = alt.Chart(states).mark_geoshape(
    fill='lightgray',
    stroke='white'
).project('albersUsa').properties(
    width=500,
    height=300
)
usa_states
Out[3]:
In [4]:
def stem_to_lat_lon(s: str) -> tuple[float, float]:
    return tuple(float(x) for x in s.split('_')[1:3])

df = pd.DataFrame([stem_to_lat_lon(fp.stem) for fp in actual], columns=['lat', 'lon'])

usa_states + alt.Chart(df).mark_circle().encode(
    longitude='lon:Q',
    latitude='lat:Q',
    size=alt.value(10),
).project(
    'albersUsa'
).properties(width=500, height=300)
Out[4]:
In [5]:
raw = pd.read_csv(folder / 'Actual_36.25_-102.95_2006_UPV_84MW_5_Min.csv')
raw['local_time'] = raw['LocalTime'].apply(lambda s: datetime.strptime(s, '%m/%d/%y %H:%M'))
raw['power'] = raw['Power(MW)']
raw = raw.drop(columns=['LocalTime', 'Power(MW)'])
raw.head()
Out[5]:
local_time power
0 2006-01-01 00:00:00 0.0
1 2006-01-01 00:05:00 0.0
2 2006-01-01 00:10:00 0.0
3 2006-01-01 00:15:00 0.0
4 2006-01-01 00:20:00 0.0
In [6]:
raw.describe()
Out[6]:
local_time power
count 105120 105120.000000
mean 2006-07-02 11:57:29.999999744 16.039188
min 2006-01-01 00:00:00 0.000000
25% 2006-04-02 05:58:45 0.000000
50% 2006-07-02 11:57:30 0.000000
75% 2006-10-01 17:56:15 32.600000
max 2006-12-31 23:55:00 84.000000
std NaN 22.735229
In [7]:
jan_1st = raw.loc[(raw['local_time'].dt.day == 1) & (raw['local_time'].dt.month == 1)]
alt.Chart(jan_1st).mark_line().encode(x='local_time', y='power', tooltip='local_time').properties(width=500)
Out[7]:
In [8]:
from datetime import date
day_avg = raw.groupby(raw['local_time'].dt.time).mean().drop(columns=['local_time']).reset_index()
day_avg['local_time'] = day_avg['local_time'].apply(lambda t: datetime.combine(date.today(), t))
alt.Chart(day_avg).mark_line().encode(x='local_time', y='power').properties(title="Average power generation over a day")
Out[8]:
In [9]:
from datetime import time

year_avg = raw.groupby(raw['local_time'].dt.date).mean().drop(columns=['local_time']).reset_index()
year_avg['local_time'] = year_avg['local_time'].apply(lambda t: datetime.combine(t, time(12, 0, 0)))
alt.Chart(year_avg).mark_line().encode(
    x='local_time', y='power'
).properties(height=500, width=800)
Out[9]:
In [10]:
sol = raw.copy()
sol['power'] /= 84  # 84 MW array, standardize on 1MW
sol['power'].mean(), len(sol)
Out[10]:
(np.float64(0.1909427094658259), 105120)
In [11]:
sol.head()
Out[11]:
local_time power
0 2006-01-01 00:00:00 0.0
1 2006-01-01 00:05:00 0.0
2 2006-01-01 00:10:00 0.0
3 2006-01-01 00:15:00 0.0
4 2006-01-01 00:20:00 0.0
In [12]:
alt.Chart(sol).mark_line().encode(
    x='local_time', y='power'
).properties(
    width=800, height=600
)
Out[12]:
In [13]:
days = raw['local_time'].dt.date.unique()[::10]
df = raw.loc[raw['local_time'].dt.date.isin(days)].assign(
    day=lambda x: x['local_time'].dt.date.apply(lambda y: datetime.combine(y, time(12, 0, 0))),
    time=lambda x: x['local_time'].dt.time.apply(lambda y: datetime.combine(date.today(), y)),
)
alt.Chart(df.sample(5000)).mark_line().encode(
    x="time:T",
    y="power:Q",
    color="day:N"
)
Out[13]:

Time to implement the uptime function. We baseline a 1MW array, then set up numerical array with loads of different sizes and batteries of different sizes. If the battery is empty, load is off. If battery is full, no chargning can occur. We measure everything in 5 minute intervals (according to the data), and assume the battery starts full. Battery state is measure in MWh stored, so in each interval we have to divide by 12.

Let's start with a naive variant (no vectorization)

In [14]:
24*365
Out[14]:
8760
In [15]:
# capacity in MWh, load in MW, sol in MW
@njit
def uptime(capacity: float, load: float, sol: list[float]) -> tuple[float, float, float, float]:
    battery: list[float] = [capacity] + [0.0 for _ in sol] # MWh
    utilization: list[float] = [0.0 for _ in sol]  # percentage
    t_interval = 24.0 * 365.0 / len(sol)  # hours
    for i in range(len(sol)):
        sol_interval: float = sol[i]
        remaining_load: float = (load - sol_interval)*t_interval
        if sol_interval > load:  # More sun than load
            utilization[i] = 1  # Can run full load
            excess_solar = sol_interval - load
            # if battery[i] < capacity:
                # battery[i+1] = battery[i] + t_interval*excess_solar
            battery[i+1] = min(battery[i] + t_interval*excess_solar, capacity)
            
        elif battery[i] > remaining_load: # Battery is full enough
            utilization[i] = 1
            battery[i+1] = battery[i] - remaining_load
        else: # Battery fully drained, cannot use all load
            utilization[i] = (sol_interval * t_interval + battery[i]) / (load * t_interval)  # equivalent: sol_interval / load + battery / (load * ts)
            battery[i+1] = 0

    batsum = 0
    for b in battery[:-1]:
        batsum += 1 if b > 0 else 0
        
    return (
        capacity,
        load,
        batsum / (len(battery) - 1),
        t_interval * sum(utilization) / (24 * 365),
    )

uptime(2.0, 1.0, sol['power'].tolist())
Out[15]:
(2.0, 1.0, 0.00023782343987823439, 0.19117101996811398)
In [16]:
@njit
def all_in_system_cost(
    solar_cost: float, battery_cost: float, load_cost: float, battery_size: float, array_size: float, sol: list[float]
) -> float:
    capacity, load, battery_util, total_util = uptime(
        capacity=battery_size / array_size,
        load=1 / array_size,
        sol=sol,
    )
    return (capacity * battery_cost + solar_cost + load_cost * load) / (load * total_util)

all_in_system_cost(200_000, 200_000, 5_000_000, 10, 5, sol['power'].tolist())
Out[16]:
10527576.740406122

Casey's number here is 1.0523e7, so I'm going to call that close enough. There appears to be a bug in Casey's implementation. The check if the battery can be charged happens at the start of the interval. If the battery is almost full this will be negative, but then there might be too much solar excess, filling the battery beyond max. However, if we leave the bug in the program, the number doesn't match quite as closely... So either this bug is not in Casey's program, or my program has another discrepancy.

In [17]:
@njit
def cost_and_elasticity(
    solar_cost: float, battery_cost: float, load_cost: float, battery_size: float, array_size: float, sol: list[float]
) -> tuple[float, float, float, float, float]:
    size_to_cost = lambda b, a: all_in_system_cost(
        solar_cost, battery_cost, load_cost, b, a, sol
    )
    cost = size_to_cost(battery_size, array_size)
    cost_battery = size_to_cost(1.01*battery_size+0.01, array_size)
    cost_array = size_to_cost(battery_size, 1.01*array_size+0.01)
    return cost, cost_battery, cost_array, (cost - cost_battery) / cost, (cost - cost_array) / cost

cost_and_elasticity(200e3, 200e3, 5e6, 10, 5, sol['power'].tolist())
Out[17]:
(10527576.740406122,
 10506720.226831798,
 10520210.612583503,
 0.0019811314691513676,
 0.0006996983260494161)

I'm running the gradient descent for more steps with a much smaller amplitude since this code appears to be ~1000 times faster than the equivalent Mathematica. This takes care of some numerical infelicities.

In [18]:
@njit
def find_minimum_system_cost(
    solar_cost: float,
    battery_cost: float,
    load_cost: float,
    sol: list[float],
) -> tuple[tuple[float, float, float], tuple[float, float, 1], tuple[float, float, float], tuple[float, float, float], any]:
    bi = min(10, 10 * load_cost / 5e6)
    ai = min(10, 10 * load_cost / 5e6)
    amplitude = 10 + 70 * (load_cost / 5e6)
    if 7e5 < load_cost < 13e5: amplitude *= 3
    if 80e6 < load_cost: amplitude *= 0.5

    cost_min = 10**10
    bi_min = bi
    ai_min = ai
    for i in range(100):
        # There was a bug here: cost_and_elasticity were getting called with wrong argument order
        # Changing them around doesn't affect the results since we always passed in the same number here
        cost, cost_bat, cost_arr, dcost_bat, dcost_arr = cost_and_elasticity(solar_cost, battery_cost, load_cost, bi, ai, sol)
        if cost < cost_min:
            ai_min, bi_min, cost_min = ai, bi, cost
        # if True: print((cost, cost_bat, cost_arr, dcost_bat, dcost_arr), bi, ai)
        bi = max(0, bi + amplitude * uniform(0.1, 1) * dcost_bat)
        ai = max(0.01, ai + amplitude * uniform(0.1, 1) * dcost_arr)

    ut = uptime(bi_min / ai_min, 1 / ai_min, sol)
    
    array_cost = solar_cost * ai_min
    storage_cost = battery_cost * bi_min
    total_cost = array_cost + storage_cost + load_cost
    
    return (
        (solar_cost, battery_cost, load_cost),
        (ai_min, bi_min, 1),
        (array_cost, storage_cost, load_cost),
        (array_cost + storage_cost, total_cost, total_cost / ut[-1]),
        ut
    )

find_minimum_system_cost(2e5, 2e5, 1e5, sol['power'].tolist())
Out[18]:
((200000.0, 200000.0, 100000.0),
 (1.3960843186370024, 0.0, 1),
 (279216.8637274005, 0.0, 100000.0),
 (279216.8637274005, 379216.8637274005, 1448531.3901021087),
 (0.0, 0.7162891142393895, 0.0, 0.26179402553414394))

That's really not that far from Casey's outcomes. I'm going to chalk the difference up to order of operations and floating point shenanigans.

In [19]:
from tqdm import tqdm
results_raw = [
    find_minimum_system_cost(200e3, 200e3, 10e3*10**(0.1*i), sol['power'].tolist())
    for i in tqdm(range(40))
]

100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 40/40 [00:16<00:00,  2.38it/s]

I'm not sure what's taking Casey's implementation so long? Let's call this a minute, then it's a factor 100? Factor 1000 if I enable compilation with numba.

In [20]:
results_raw
Out[20]:
[((200000.0, 200000.0, 10000.0),
  (1.2244353216573673, 0.0, 1),
  (244887.06433147346, 0.0, 10000.0),
  (244887.06433147346, 254887.06433147346, 1091908.353259966),
  (0.0, 0.8167029995887599, 0.0, 0.23343265354687595)),
 ((200000.0, 200000.0, 12589.254117941673),
  (1.2407553998584486, 0.0, 1),
  (248151.07997168973, 0.0, 12589.254117941673),
  (248151.07997168973, 260740.3340896314, 1102931.8827598882),
  (0.0, 0.8059606269810189, 0.0, 0.2364065616066658)),
 ((200000.0, 200000.0, 15848.931924611135),
  (1.2426763817356783, 0.0, 1),
  (248535.27634713566, 0.0, 15848.931924611135),
  (248535.27634713566, 264384.2082717468, 1116702.2814634265),
  (0.0, 0.8047147388471921, 0.0, 0.23675442654713136)),
 ((200000.0, 200000.0, 19952.623149688796),
  (1.2574119765816938, 0.0, 1),
  (251482.39531633875, 0.0, 19952.623149688796),
  (251482.39531633875, 271435.01846602757, 1133788.5699939742),
  (0.0, 0.7952842971311004, 0.0, 0.2394053226938689)),
 ((200000.0, 200000.0, 25118.8643150958),
  (1.2745379278457512, 0.0, 1),
  (254907.58556915022, 0.0, 25118.8643150958),
  (254907.58556915022, 280026.44988424605, 1155045.7807352352),
  (0.0, 0.7845980713106118, 0.0, 0.2424375332603678)),
 ((200000.0, 200000.0, 31622.776601683796),
  (1.2907676411166946, 0.0, 1),
  (258153.52822333892, 0.0, 31622.776601683796),
  (258153.52822333892, 289776.3048250227, 1181508.0168288313),
  (0.0, 0.7747327777250909, 0.0, 0.24525970259836458)),
 ((200000.0, 200000.0, 39810.717055349734),
  (1.3074071903686766, 0.0, 1),
  (261481.4380737353, 0.0, 39810.717055349734),
  (261481.4380737353, 301292.15512908506, 1214449.519754507),
  (0.0, 0.764872648220643, 0.0, 0.24808948435336303)),
 ((200000.0, 200000.0, 50118.723362727236),
  (1.3242653271693425, 0.0, 1),
  (264853.0654338685, 0.0, 50118.723362727236),
  (264853.0654338685, 314971.78879659576, 1255455.1940033447),
  (0.0, 0.755135681259231, 0.0, 0.25088254069205485)),
 ((200000.0, 200000.0, 63095.734448019335),
  (1.3441920069020703, 0.0, 1),
  (268838.4013804141, 0.0, 63095.734448019335),
  (268838.4013804141, 331934.13582843344, 1306479.9748157333),
  (0.0, 0.7439413378931466, 0.0, 0.2540675266570769)),
 ((200000.0, 200000.0, 79432.82347242816),
  (1.3679954495180282, 0.0, 1),
  (273599.08990360564, 0.0, 79432.82347242816),
  (273599.08990360564, 353031.9133760338, 1369879.5997687874),
  (0.0, 0.7309965836160638, 0.0, 0.25771017645318584)),
 ((200000.0, 200000.0, 100000.0),
  (1.3960842407901883, 0.0, 1),
  (279216.84815803764, 0.0, 100000.0),
  (279216.84815803764, 379216.84815803764, 1448531.39152861),
  (0.0, 0.716289154180264, 0.0, 0.26179401452795353)),
 ((200000.0, 200000.0, 125892.54117941674),
  (1.4272348255765002, 0.0, 1),
  (285446.96511530003, 0.0, 125892.54117941674),
  (285446.96511530003, 411339.5062947168, 1546021.2933650373),
  (0.0, 0.7006555488134701, 0.0, 0.2660632864890262)),
 ((200000.0, 200000.0, 158489.3192461114),
  (1.462893346767889, 0.0, 1),
  (292578.66935357783, 0.0, 158489.3192461114),
  (292578.66935357783, 451067.9885996892, 1666774.9112474346),
  (0.0, 0.683576832316174, 0.0, 0.27062321706180736)),
 ((200000.0, 200000.0, 199526.23149688798),
  (1.5049548149112784, 0.0, 1),
  (300990.9629822557, 0.0, 199526.23149688798),
  (300990.9629822557, 500517.19447914365, 1816227.4003583568),
  (0.0, 0.664471776887835, 0.0, 0.2755806868569362)),
 ((200000.0, 200000.0, 251188.6431509581),
  (1.5581977132781593, 0.0, 1),
  (311639.54265563184, 0.0, 251188.6431509581),
  (311639.54265563184, 562828.1858065899, 2000871.4452022547),
  (0.0, 0.6417670822377125, 0.0, 0.28129152782711503)),
 ((200000.0, 200000.0, 316227.7660168379),
  (1.6374665699018067, 0.0, 1),
  (327493.31398036133, 0.0, 316227.7660168379),
  (327493.31398036133, 643721.0799971992, 2228338.2683748086),
  (0.0, 0.6106994905306473, 0.0, 0.28887942604274514)),
 ((200000.0, 200000.0, 398107.17055349733),
  (1.717843162614406, 0.0, 1),
  (343568.6325228812, 0.0, 398107.17055349733),
  (343568.6325228812, 741675.8030763785, 2508317.714592094),
  (0.0, 0.5821253195653135, 0.0, 0.29568654670885297)),
 ((200000.0, 200000.0, 501187.23362727254),
  (1.8188832686304455, 0.001872700053609922, 1),
  (363776.6537260891, 374.5400107219844, 501187.23362727254),
  (364151.1937368111, 865338.4273640837, 2852194.254691209),
  (0.0010295878168256496,
   0.549787893069666,
   0.17863394216133943,
   0.30339393116047386)),
 ((200000.0, 200000.0, 630957.3444801933),
  (1.964273686206744, 0.020142225802007474, 1),
  (392854.73724134883, 4028.445160401495, 630957.3444801933),
  (396883.18240175035, 1027840.5268819437, 3273285.8998570726),
  (0.010254286835611288,
   0.5090940264699692,
   0.20901826484018265,
   0.31400878454485875)),
 ((200000.0, 200000.0, 794328.2347242817),
  (2.0767665656975094, 0.03176161149196756, 1),
  (415353.3131395019, 6352.322298393512, 794328.2347242817),
  (421705.6354378954, 1216033.8701621771, 3788682.124269707),
  (0.015293780252717042,
   0.48151776733950685,
   0.22544710806697107,
   0.32096487123383977)),
 ((200000.0, 200000.0, 1000000.0),
  (2.2657804276118165, 0.13158885393762265, 1),
  (453156.0855223633, 26317.77078752453, 1000000.0),
  (479473.85630988784, 1479473.856309888, 4420033.582871287),
  (0.05807661339731858,
   0.44134903268364023,
   0.2583618721461187,
   0.33472004874424743)),
 ((200000.0, 200000.0, 1258925.4117941675),
  (2.6786303117184183, 1.2519415977284671, 1),
  (535726.0623436837, 250388.31954569343, 1258925.4117941675),
  (786114.3818893771, 2045039.7936835447, 5206619.439035858),
  (0.46738125535707487,
   0.37332512651156824,
   0.34356925418569256,
   0.39277689057724513)),
 ((200000.0, 200000.0, 1584893.1924611141),
  (3.2072796988205106, 3.254418751601248, 1),
  (641455.9397641021, 650883.7503202496, 1584893.1924611141),
  (1292339.6900843517, 2877232.8825454656, 5985962.389398696),
  (1.0146975185226512,
   0.31179070549031124,
   0.4385369101978691,
   0.48066337463815745)),
 ((200000.0, 200000.0, 1995262.3149688807),
  (3.724585678506174, 4.935711114499921, 1),
  (744917.1357012348, 987142.2228999842, 1995262.3149688807),
  (1732059.358601219, 3727321.6735700998, 6726135.534941762),
  (1.3251705130540843,
   0.26848623882403794,
   0.5173611111111112,
   0.5541550053826522)),
 ((200000.0, 200000.0, 2511886.4315095823),
  (4.368349264325126, 7.289540329291121, 1),
  (873669.8528650252, 1457908.0658582242, 2511886.4315095823),
  (2331577.9187232493, 4843464.350232832, 7427685.756192738),
  (1.6687173777111652,
   0.228919424590582,
   0.6210140791476408,
   0.6520825610042341)),
 ((200000.0, 200000.0, 3162277.6601683795),
  (5.21958033649443, 10.534271785180508, 1),
  (1043916.067298886, 2106854.3570361016, 3162277.6601683795),
  (3150770.4243349875, 6313048.084503368, 8056370.036931467),
  (2.0182219845389957,
   0.1915862838642731,
   0.7587899543378995,
   0.7836094985165178)),
 ((200000.0, 200000.0, 3981071.705534973),
  (6.099504721509909, 13.246110297384787, 1),
  (1219900.944301982, 2649222.059476957, 3981071.705534973),
  (3869123.003778939, 7850194.709313912, 8791140.433669169),
  (2.1716698161856267,
   0.16394773767015855,
   0.8775399543378996,
   0.8929665915981128)),
 ((200000.0, 200000.0, 5011872.336272725),
  (6.706569334512411, 13.634821884925795, 1),
  (1341313.8669024822, 2726964.3769851592, 5011872.336272725),
  (4068278.2438876415, 9080150.580160366, 9931826.565222971),
  (2.033054637154078,
   0.14910753175307376,
   0.9028443683409437,
   0.9142477992875135)),
 ((200000.0, 200000.0, 6309573.444801937),
  (6.957547402495498, 13.905072283923127, 1),
  (1391509.4804990997, 2781014.4567846255, 6309573.444801937),
  (4172523.9372837255, 10482097.382085662, 11345369.224256417),
  (1.9985594749862179,
   0.14372880875253902,
   0.914630898021309,
   0.9239097622027956)),
 ((200000.0, 200000.0, 7943282.347242822),
  (7.2203969015187415, 14.250950156232228, 1),
  (1444079.3803037482, 2850190.0312464456, 7943282.347242822),
  (4294269.411550194, 12237551.758793015, 13102838.429597715),
  (1.9737073114685257,
   0.1384965416222008,
   0.9263698630136986,
   0.93396189112352)),
 ((200000.0, 200000.0, 10000000.0),
  (7.526708924007613, 14.545195378002349, 1),
  (1505341.7848015225, 2909039.07560047, 10000000.0),
  (4414380.860401993, 14414380.860401992, 15291825.076234654),
  (1.9324774645673062,
   0.1328601929603447,
   0.9364345509893455,
   0.9426200462365792)),
 ((200000.0, 200000.0, 12589254.117941676),
  (7.977919068714, 14.781168116896831, 1),
  (1595583.8137428, 2956233.623379366, 12589254.117941676),
  (4551817.437122166, 17141071.555063844, 18025501.781149704),
  (1.8527598474722908,
   0.1253459694673482,
   0.9459094368340943,
   0.9509345017506942)),
 ((200000.0, 200000.0, 15848931.92461114),
  (8.268734937956264, 14.986826159649816, 1),
  (1653746.987591253, 2997365.231929963, 15848931.92461114),
  (4651112.219521216, 20500044.144132357, 21440885.925713792),
  (1.8124690502358785,
   0.12093748408957516,
   0.9518930745814308,
   0.9561192674201444)),
 ((200000.0, 200000.0, 19952623.14968881),
  (8.45676340044536, 15.108700439348954, 1),
  (1691352.6800890719, 3021740.087869791, 19952623.14968881),
  (4713092.767958863, 24665715.917647675, 25722453.371079743),
  (1.7865819018364975,
   0.11824854884166879,
   0.9552321156773211,
   0.958917703603647)),
 ((200000.0, 200000.0, 25118864.31509582),
  (9.399497527776898, 15.205708971058254, 1),
  (1879899.5055553797, 3041141.7942116507, 25118864.31509582),
  (4921041.29976703, 30039905.61486285, 31087424.98987767),
  (1.6177150880803093,
   0.10638866567546328,
   0.9639840182648401,
   0.966304080336152)),
 ((200000.0, 200000.0, 31622776.601683795),
  (10.261365548131185, 15.311778322952813, 1),
  (2052273.1096262368, 3062355.6645905627, 31622776.601683795),
  (5114628.774216799, 36737405.3759006, 37794221.796346754),
  (1.492177454465738,
   0.09745291650604158,
   0.9702054794520548,
   0.9720376192387082)),
 ((200000.0, 200000.0, 39810717.05534973),
  (11.01483739945717, 15.395388021172348, 1),
  (2202967.479891434, 3079077.6042344696, 39810717.05534973),
  (5282045.084125903, 45092762.139475636, 46190757.520774364),
  (1.3976954414171434,
   0.09078663295104854,
   0.9747336377473363,
   0.976229110752191)),
 ((200000.0, 200000.0, 50118723.36272725),
  (12.088155833176115, 15.526193350706933, 1),
  (2417631.166635223, 3105238.670141387, 50118723.36272725),
  (5522869.83677661, 55641593.199503854, 56712563.20921373),
  (1.2844137323325262,
   0.08272560461666831,
   0.9799752663622526,
   0.981115824270558)),
 ((200000.0, 200000.0, 63095734.44801936),
  (13.26337283286767, 15.631386219535369, 1),
  (2652674.566573534, 3126277.2439070735, 63095734.44801936),
  (5778951.810480608, 68874686.25849997, 69889228.88845776),
  (1.178537798530369,
   0.07539560356185737,
   0.9844939117199392,
   0.9854835624016257)),
 ((200000.0, 200000.0, 79432823.47242822),
  (14.579004849782326, 15.83435752110239, 1),
  (2915800.969956465, 3166871.5042204782, 79432823.47242822),
  (6082672.474176943, 85515495.94660516, 86419058.84507793),
  (1.0861068834433378,
   0.06859178732044462,
   0.988898401826484,
   0.9895444024669074))]
In [21]:
labels = [
    'solar cost ($/MW)', 'battery cost ($/MW)', 'load cost ($/MW)',
    'array size (MW)', 'battery size (MWh)', 'load size (1 MW by definition)',
    'array cost ($)', 'battery cost ($)', 'load cost ($, normalized to 1 MW)',
    'total power system cost ($)', 'total system cost ($)', 'total cost per utilization ($)',
    'battery size relative to 1 MW array', 'load size relative to 1 MW array', 'annual battery utilization', 'annual load utilization'
]

results = pd.DataFrame([dict(zip(labels, (x for tup in r for x in tup))) for r in results_raw])
results.head()
Out[21]:
solar cost ($/MW) battery cost ($/MW) load cost ($/MW) array size (MW) battery size (MWh) load size (1 MW by definition) array cost ($) battery cost ($) load cost ($, normalized to 1 MW) total power system cost ($) total system cost ($) total cost per utilization ($) battery size relative to 1 MW array load size relative to 1 MW array annual battery utilization annual load utilization
0 200000.0 200000.0 10000.000000 1.224435 0.0 1 244887.064331 0.0 10000.000000 244887.064331 254887.064331 1.091908e+06 0.0 0.816703 0.0 0.233433
1 200000.0 200000.0 12589.254118 1.240755 0.0 1 248151.079972 0.0 12589.254118 248151.079972 260740.334090 1.102932e+06 0.0 0.805961 0.0 0.236407
2 200000.0 200000.0 15848.931925 1.242676 0.0 1 248535.276347 0.0 15848.931925 248535.276347 264384.208272 1.116702e+06 0.0 0.804715 0.0 0.236754
3 200000.0 200000.0 19952.623150 1.257412 0.0 1 251482.395316 0.0 19952.623150 251482.395316 271435.018466 1.133789e+06 0.0 0.795284 0.0 0.239405
4 200000.0 200000.0 25118.864315 1.274538 0.0 1 254907.585569 0.0 25118.864315 254907.585569 280026.449884 1.155046e+06 0.0 0.784598 0.0 0.242438
In [22]:
subsystems = results.melt(id_vars='annual load utilization', value_vars=[
    'array cost ($)', 'battery cost ($)', 'load cost ($/MW)', 'total power system cost ($)', 'total system cost ($)', 
])
subsystems = subsystems.loc[subsystems['value'] <= 2e7]
alt.Chart(subsystems).mark_line().encode(
    x='annual load utilization:Q',
    y=alt.Y('value:Q', scale=alt.Scale(domain=[0, 2e7], clamp=True)),
    color='variable:N',
    tooltip='variable:N',
).properties(width=800, height=600)
Out[22]:
In [23]:
subsystems = results.melt(id_vars="load cost ($/MW)", value_vars=[
    'array cost ($)', 'battery cost ($)', 'load cost ($/MW)', 'total power system cost ($)', 'total system cost ($)', 
    'total cost per utilization ($)',
])
subsystems = subsystems.loc[subsystems['value'] > 0]

alt.Chart(subsystems).mark_line().encode(
    x=alt.X('load cost ($/MW):Q', scale=alt.Scale(type='log', domain=[0.5e4, 1e8])),
    y=alt.Y('value:Q', scale=alt.Scale(type='log', domain=[0.5e3, 1e8])),
    color='variable:N',
    tooltip=('variable:N','value:Q'),
).properties(width=800, height=600)
Out[23]:
In [24]:
alt.Chart(results).mark_line().encode(
    x=alt.X('load cost ($/MW)', scale=alt.Scale(type='log')),
    y='annual load utilization'
).properties(width=800, height=600)
Out[24]:
In [25]:
results.describe()
Out[25]:
solar cost ($/MW) battery cost ($/MW) load cost ($/MW) array size (MW) battery size (MWh) load size (1 MW by definition) array cost ($) battery cost ($) load cost ($, normalized to 1 MW) total power system cost ($) total system cost ($) total cost per utilization ($) battery size relative to 1 MW array load size relative to 1 MW array annual battery utilization annual load utilization
count 40.0 40.0 4.000000e+01 40.000000 40.000000 40.0 4.000000e+01 4.000000e+01 4.000000e+01 4.000000e+01 4.000000e+01 4.000000e+01 40.000000 40.000000 40.000000 40.000000
mean 200000.0 200000.0 9.654325e+06 4.518091 5.870373 1.0 9.036181e+05 1.174075e+06 9.654325e+06 2.077693e+06 1.173202e+07 1.311563e+07 0.754919 0.433660 0.420597 0.546628
std 0.0 0.0 1.851465e+07 3.932988 6.983019 0.0 7.865975e+05 1.396604e+06 1.851465e+07 2.151899e+06 2.023223e+07 2.006034e+07 0.854375 0.284317 0.431302 0.322336
min 200000.0 200000.0 1.000000e+04 1.224435 0.000000 1.0 2.448871e+05 0.000000e+00 1.000000e+04 2.448871e+05 2.548871e+05 1.091908e+06 0.000000 0.068592 0.000000 0.233433
25% 200000.0 200000.0 9.485821e+04 1.389062 0.000000 1.0 2.778124e+05 0.000000e+00 9.485821e+04 2.778124e+05 3.726706e+05 1.428868e+06 0.000000 0.137087 0.000000 0.260773
50% 200000.0 200000.0 8.971641e+05 2.171273 0.081675 1.0 4.342547e+05 1.633505e+04 8.971641e+05 4.505897e+05 1.347754e+06 4.104358e+06 0.036685 0.461433 0.241904 0.327842
75% 200000.0 200000.0 8.457462e+06 7.296975 14.324511 1.0 1.459395e+06 2.864902e+06 8.457462e+06 4.324297e+06 1.278176e+07 1.365009e+07 1.630466 0.719966 0.928886 0.936126
max 200000.0 200000.0 7.943282e+07 14.579005 15.834358 1.0 2.915801e+06 3.166872e+06 7.943282e+07 6.082672e+06 8.551550e+07 8.641906e+07 2.171670 0.816703 0.988898 0.989544
In [26]:
for col in ['array cost ($)', 'battery cost ($)', 'total power system cost ($)']:
    results[col.replace('$', '$/MWh')] = results[col] / (10 * 24 * 365 * results['annual load utilization'])

import numpy as np
min_util = results.loc[results['annual load utilization'].idxmin()]
uvals = np.arange(0.01, min_util['annual load utilization'], 0.001)
underutilized_solar = pd.Series(uvals, name='annual load utilization').to_frame()
underutilized_solar['variable'] = 'underutilized solar'
underutilized_solar['value'] = min_util['array cost ($)'] / (10 * 24 * 365 * uvals)

subsystems = results.melt(id_vars='annual load utilization', value_vars=[
    'array cost ($/MWh)', 'battery cost ($/MWh)', 'total power system cost ($/MWh)', 
])
subsystems = pd.concat([subsystems, underutilized_solar])
alt.Chart(subsystems).mark_line().encode(
    x='annual load utilization:Q',
    y=alt.Y('value:Q', scale=alt.Scale(domain=[0, 80])),
    color='variable:N',
).properties(width=800, height=600).interactive()
Out[26]:
In [ ]: